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ABSTRACT 

We describe a made-to-measure algorithm for constructing TV-particle models of stel- 
lar systems from observational data (x 2 M2M), extending earlier ideas by Syer and 
Trcmaine. The algorithm properly accounts for observational errors, is flexible, and 
can be applied to various systems and geometries. We implement this algorithm in a 
parallel code NMAGIC and carry out a sequence of tests to illustrate its power and 
performance: (i) We reconstruct an isotropic Hernquist model from density moments 
and projected kinematics and recover the correct differential energy distribution and 
intrinsic kinematics, (ii) We build a self-consistent oblate three- integral maximum ro- 
tator model and compare how the distribution function is recovered from integral held 
and slit kinematic data, (iii) We create a non-rotating and a figure rotating triaxial 
stellar particle model, reproduce the projected kinematics of the figure rotating system 
by a non-rotating system of the same intrinsic shape, and illustrate the signature of 
pattern rotation in this model. From these tests we comment on the dependence of the 
results from x 2 M2M on the initial model, the geometry, and the amount of available 
data. 
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1 INTRODUCTION 

Understanding the structure and dynamics of galaxies re- 
quires knowledge of the total gravitational potential and 
the distribution of stellar orbits. Due to projection effects 
the orbital structure is not directly given by observations. 
In equilibrium stellar systems, the phase-space distribution 
function (DF) fully determines the state of the galaxy. Dy- 
namical models of observed galaxies attempt to recover their 
DF and total (i.e. due to visible and dark matter) grav- 
itational potential consistent with the observational data. 
Several m ethods to tackle this prob lem exist. Jean's theo- 
rem (e.g. iBinnev fc Tremaind Il987r i requires that the DF 
depends on the phase-space coordinates only through the 
integrals of motion. If these integrals can be expressed 
or approximated in terms of analytic functions, one can 
parametrize the DF explicitly. This approach ha s been ap- 
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physical reason why the DF should only depend on the clas- 
sical integrals and most orbits in axisymmetric systems have 
an approxi mate third integr al of motion, which is not known 
in general dOuongrenl l l962l ). 

ISchwarzschildl (|l979h developed a technique for nu- 
merically building self-consistent models of galaxies, with- 
out explicit knowledge of the integrals of motion. In 
this method, a library of orbits is computed and or- 
bits are then superposed with positive definite weights 
to reproduce observed photometry and kinematics. The 
Schwarzschild method has been used to model stellar sys- 
tems for measurements of global mass-to-light ratios, in- 
ternal kinematics an d the masses of central supermas- 
sive black holes (e.g. iRix et all [l997l : ICretton et"ahl Il999l : 
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Romanowskv fc Kochanek 2003; ICappellari et al. 20021: 



Verolme et al.ll2002l; iGebhardt et al.ll2003l; Ivan de Ven et al 



2005 : IValluri et alJ I2004L ICopin et al J 12001 iThomas et al 
20051 ). The method is well-tested, and modern implemen- 
tations are quite efficient. However, it also has some draw- 
backs: symmetry assumptions are often made, and the po- 
tential must be chosen a priori. Initial conditions for a rep- 
resentative orbit library have to be carefully chosen, which 
becomes more complicated as the complexity of the poten- 
tial's phase space structure increases, in terms of number of 
orbit families, resonances, chaotic and semi-chaotic regions. 
As a result, most Schwarzschild models in the literature to 
date are axisymmetric. 

Thus there is scope for e xploring alternative ap- 
proaches. ISver fc Tremainel (1 19961 . hereafter ST96) invented 
a particle-based algorithm for constructing models of stellar 
systems. This "made-to-measure" (M2M) method works by 
adjusting individually adaptable weights of the particles as 
a function of time, until the model converges to the obser- 
vational data. The first practical application of the M2M 
method constructed a d ynamical model of th e Milky Way's 
barred bulge and disk (Biss antz et al.l 2004) and was able 
to match the event timescale distribution of microlensing 
events towards the bulge. This paper illustrates some of the 
promise that lies in particle-based methods, in that it was 
relatively easy to model a rapidly rotating stellar system. 
However, other important modeling aspects were not yet 
implemented, such as a proper treatment of observational 
errors. The purpose of the present paper is to show how this 
can be done, and to describe and test our modified x 2 M2M 
method designed for this purpose. 

The paper is organized as follows. In the Section [2] we 
describe the M2M algorithm of ST96. Then in Section[3]we 
extend the algorithm in order to include observational er- 
rors. We also discuss how we include density and kinematic 
observables in the same model, and describe the NMAGIC 
code, our parallel implementation of the x 2 M2M method. In 
Section [4] we present the models we use to test this imple- 
mentation, and the results of these tests follow in Section [S] 
Finally, the paper closes with the conclusions in Section [6] 



2 SYER & TREMAINE'S 

MADE-TO-MEASURE ALGORITHM 

The M2M algorithm is designed to build a particle model to 
match the observables of some target system. The algorithm 
works by varying the individually adaptable weights of the 
particles moving in the global potential until the model min- 
imizes deviations between its observables and those of the 
target. An observable of a system characterized by a distri- 
bution function /(z), is defined as 



Y 3 = / ^(z)/(z) d 6 



(1) 



where Kj is a known kernel and z = (r, v) are phase-space 
coordinates. Examples of typical observables include surface 
or volume densities and line-of-sight kinematics. The equiv- 
alent observable of the particle model is given by 



JV 



Vi(t) = ^WiKj [z, 



it)]: 



(2) 



where it), are the weights and Zi are the phase-space coordi- 
nates of the particles, i = 1, • • • , N. In the following, we use 
units and normalization such that 



(3) 



so that the equivalent masses of the particles are TOj = WiM 
with M the total mass of the system. 

Given a set of observables Yj, j = 1, • ■ • , J, we want to 
construct a system of N particles i — 1, • • • , N orbiting in 
the potential, such that the observables of the system match 
those of the target system. The heart of the algorithm is a 
prescription for changing particle weights by specifying the 
"force-of-change" (hereafter FOC): 



dt 



Here 



Aj(t) 



3 



vAt) 



(4) 



(5) 



measures the deviation between target and model observ- 
ables. The constant e is small and positive and, to this point, 
the Zj are arbitrary constants. The linear dependence of 
the FOC for weight Wi on Wi itself ensures that the particle 
weights cannot become negative, and the dependence on the 
kernel Kj ensures that a mismatch in observable j only has 
influence of the weight of particle i when that particle actu- 
ally contributes to the observable j. The choice of A in terms 
of the ratio of the model and targe t obse rvables makes the 
algorithm closely related to Lucy's l|l974h method, in which 
one iteratively solves an integral equation for the distribu- 
tion underlying the process from observational data. 

Since in typical applications the number of particles 
greatly exceeds the number of independent constraints, the 
solutions of the set of differential equations Q are under- 
determined, i.e. the observables of the particle model can 
remain constant, even as the particle weights may still be 
changing with time. To remove this ill-conditioning, ST96 
maximize the function 



with 



and the entropy 

S = — Wi \og(wi/wi) 



(6) 



(7) 



(8) 



as a profit function. The {u>i} are a predetermined set of 
weights, the so-called priors. Since 

dS 

— = -M(l°g(w«M) + !). ( 9 ) 
dun 

if a particle weight Wi < Wife then equation (|9]) becomes 
positive while it is negative when Wi > ibi/e. Therefore the 
entropy term pushes the particle weights to remain close to 
their priors (more specifically, close to Wi/e). Equation Q 
is now replaced by 
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E 



Kj [gift)] 



A f (t) 



(10) 



If we replace the FOC equation (|10[) by 



with Zj now fixed to Yj by the requirement that equation ((6} 
will be maximized, as discussed in Section [3] The constant 
/x governs the relative importance of the entropy term in 
equation (1101) : When \i is large the {wi} will remain close 
to their priors {u>i}. In the following, we will generally set 
Wi — wo = 1/N; i.e., the particle distribution follows the 
initial model, but this is not necessary. 

To reduce temporal fluctuations, ST96 introduced tem- 
poral smoothing by substituting Aj (t) in Equations @ and 
JT01) with 



(11) 



A j {t) = a A 3 (t-r)e- aT dr,, 
Jo 

which can be expressed in the form of the differential equa- 
tion 



(12) 



The smoothing time is 1/a. The temporal smoothing sup- 
presses fluctuations in the model observables and hence in 
the FOC correction of the particle weights - in the compu- 
tation of these quantities the effective number of particles is 
increased as each particle is effectively smeared backwards 
in time along its orbit. The smoothing time should satisfy 
2e < a to avoid excessive temporal smoothing, which slows 
down convergence. 



3 x -BASED MADE-TO-MEASURE 
ALGORITHM TO MODEL 
OBSERVATIONAL DATA 

The M2M algorithm as originally formulat ed by ST96 is 
well adapted to modeling density fields (e.g. iBissantz et al.l 
I2004T ). It is not, however, well suited to mixed observables 
such as densities and kinematics, where the various ratios of 
model to target observable can take widely different values, 
or to problems where observables can become zero, when A 
diverges. Moreover, the \ 2 defined as in equations (|7I5|) is 
not the usual one, but is given by the relative deviations 
between model and data. Thus extremizing F (equation [SJ) 
with this x 2 does not produce the best model given the 
observed data. We have therefore modified the M2M method 
as described in this section. 

We begin by considering observational errors. We do 
this by replacing equation © by 



A*(*) = 



Vj 



(13) 



where <r(Yj) in the denominator is the error in the target 
observable. With this definition of Aj equation (0 now mea- 
sures the usual absolute x 2 - As a result of this, if we now 
maximize the function of equation^ with respect to the Wi's 
we obtain the condition 



dS 
' dwi 



E- 



Kji 



(Yj. 



: 0. 



(14) 



duui(t) 
dt 



£Wi(t) fj, 



dS ^ Kj fc(t)] 

2f a{Yj) ^ 



(15) 



then the particle weights will have converged once F is max- 
imized with respect to all Wi, i.e., once the different terms 
in the bracket balance. For large fi, the solutions of eqn. [15] 
will have smooth weight distributions at the expense of a 
compromise in matching \ 2 ■ 

In the absence of the entropy term, the solutions of 
eqs. [15] near convergence can be characterized by an argu- 
ment closely similar to that used by ST96 to study the solu- 
tions of their eqs. (4). For small e, the weights Wi(t) change 
only over many orbits, so we can orbit-average over periods 
iorb < r < torb/e and write the equations for the orbit- 
averaged < Aj > as 



d<A 3 (t) > 
dt 



eAjk < Afc(t) >, 



where the matrix A has components 

< Kjj>< K ki > o 

A.jk — i-'i Wi , 

O-jCTk 



(16) 



(17) 



and we have replaced Wi(t) by the constant ui^, because near 
convergence the dominant time-dependence is in < A*, > 
rather than Wi. The matrix A is symmetric by construction 
and positive definite, i.e., x -A-x > for all vectors x; so all 
its eigenvalues are real and positive. The solutions to eqs. 1161 
then converge exponentially to < Aj(t) >= 0. As for eqs. (4) 
of ST96, this argument suggests that if e is sufficiently small 
and we start close to the correct final solution, then the 
model observables converge to their correct final values on 
0(e~ 1 ) orbital periods. 

Substituting Aj in equation pip leads to 



AAt) 



Vi(t)-Yj 



(18) 



which allows us to temporally smooth model observables 
directly 



Vj{t) = a / yj(t- r)e QT dr. 
Jo 



(19) 



This corrects the typo in equation (19) of ST96. 



In practice, yj can be computed using the equivalent differ- 
ential equation, in the same manner as before. 

Since the uncertainty in any observable presumably 
never becomes zero, the Aj in equation (|13|l remain well- 
defined even when the observables themselves take zero val- 
ues. However, if the data entering \ 2 have widely different 
relative errors, the FOC equation may be dominated by only 
a few of the Aj . This can slow down convergence of the other 
observables and thus lead to noisy final models. Also, notice 
that the cost of deriving the FOC from minimizing \ 2 is that 
equation ([6]) is maximized only if the observables are exactly 
of the form given by equation i.e. the kernel Kij may 
depend on the particle's phase-space coordinates but must 
not depend on its weight Wi. 

We adopt the convention throughout this paper in 
which the positive x-axis points in the direction of the ob- 
server, so that a particle with velocity v x < will be moving 
away from the observer. 

Our implementation of the x 2 M2M algorithm models 
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volume luminosity densities (equivalent to luminous mass 
densities for constant mass-to-light ratio), and line-of-sight 
velocities. As in the Schwarzschild method, dark matter, 
which generally has a different spatial distribution from the 
stars, can be included as an external potential, to be added 
to the potential from the luminous particles. The form of the 
dark matter potential can be guided by cosmological simu- 
lations, or also include information from gas velocities and 
other data. 



3.1 Densities 

For modeling the target distribution of stars one can use as 
observables the surface density or space density in various 
grids, or also some functional representations such as, e.g. , 
isophote fits, multi-Gaussian expansions, etc. In this paper 
we have chosen to model a spherical harmonics expansion of 
the three-dimensional density, where we expand the density 
in surface harmonics computed on a 1-D radial mesh of radii 
rfe. The expansion coefficients, Aim are computed based on 
a cloud-in-cell scheme. The function 



' r h+1 -r k 





if r € [r fc _i,r fe ) 
if r e [r k ,r k+1 ] 
otherwise, 



gives the fractional contribution of the weight w of a particle 
at radius r to shell k. The model observable is then the mass 
on each shell k, 



m k = Mj2 w ^k IC (n) = MYwaZ IC 



(20) 



Comparing with equation (2), we recognize the kernel for 
this observable as 



K ki = M>yZ Ic . 



(21) 



Thus the FOC on a particle is computed by linear interpo- 
lation of the contributions from the adjacent shells. From 
equation (|13[) . we obtain 



A k [m] = 



m k — M k 



(22) 



a{M k ) 

where M k is the target mass on shell k and a(M k ) its un- 
certainty. 

The spherical harmonic coefficients for the particle 
model with I > are computed via 



ai m ,k = M 7ki C Yi rn (6i, ipi)wj 



Now the kernel is given by 

K» = Mi& IC Yr(Bi, Vi), J = {lm, k}, 



(23) 



(24) 



and depends on the spherical harmonics; the same variation 
also holds therefore for the FOC. From equation (|13p . we 
obtain 



» r i Q>lm,k Aim .k . n 7 -i 

Aj[a !m ] = — r — , j = {lm,k\, 

cr{Ai mi k) 



(25) 



with Aim,k as the target moments and a(Ai m ,k) as their 
errors. aoo,fc and Aoo.fc are of course related to the mass on 
shell k via the relation \/47raoo.fc = rn k , etc., but we will use 
the masses on shells m k , M k as observables in the following. 



3.2 Kinematics 

Unlike for the density observables, we use kinematic observ- 
ables computed in the plane of the sky to compare with the 
target model. Since kinematic data can either come from 
restricted spatial regions {e.g. slit spectra) or from integral 
fields, we do not specify any special geometry for computing 
these observables. 

The shape of the line-of-sight velocity distribution 
(LOSVD) can be expressed in a truncated Gauss-Hermite 
series with coefficients h n , n = 1, • • • , n max (iGerhard 1 ll993l . 
Ivan der Marel fc Franxl [l993). Since the kernel in equa- 
tion (115[) cannot depend on masses, this puts some con- 
straints on which observables can be used in the FOC. 
For kinematics, suitable observables are the mass-weighted 
Gauss-Hermite coefficients, which we use as follows. Particle 
weights are assigned to a spatial cell, C p , of the kinematic 
observable under consideration using the selection function 



&pi — 



1 if (yi,Zi) £ C p 
otherwise. 



This selection function can be replaced appropriately if see- 
ing conditions need to be taken into account. In our present 
application this is not necessary. The mass-weighted kine- 
matic moments are computed as 



m p h n ,p = 2\fwM 8piU n (vpi)Wi 



Vpi — {V x ,i — V p ) j a v , 



(26) 



(27) 



and where m p is the mass in cell C p , and the dimensionless 
Gauss-Hermite functions (jGerhardl 1993T ) 



(y) = (2 n+1 7rn!) 1/2 H n (u) exp (-v 2 /2) 



(28) 



H n are the standard Hermite polynomials. For the mass- 
weighted higher order moments we obtain the kernel 



Kji = 2^ c KMSp i u n (vpi), j = {n,p}. 



and as usual 



Aj[m hn] 



a(B n! p) 



{n,p}. 



(29) 



(30) 



The velocity V p and dispersion a p are not free parameters; 
rather we set V p and a p to the mean line-of-sight veloc- 
ity and velocity dispersion obtained from the best fitting 
Gaussian to the observed (target) LOSVD. This implies 
Bi )P = (m p fei, p )targct = B2, P = (m p /i2,p)tar g ot = for the 
first and second order mass-weighted target Gauss-Hermite 
coefficients. If the model b 1>p and b2,p both converge to zero, 
then the LOSVD of the particle model automatically has the 
correct mean line-of-sight velocity and velocity dispersion. 
For describing the higher-order structure of the LOSVD we 
include terms m h n (n — 1, • • • , 4) in the test modeling de- 
scribed below. 



3.3 Implementation: the NMAGIC parallel code 

The routine for updating the particle weights includes three 
main steps: First, all the observables used in the modeling 
process are computed as described above. Then we change 
the particle weights in accordance with equation (|15[) by 
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Initial conditions 




Stop 



Figure 1. A high level flowchart describing NMAGIC. The main x 2 M2M algorithm is contained in the dashed block, the remainder is 
an optional potential solver and code for moving the particles, both of which are exchangeable. In our tests, x 2 M2M is generally applied 
only after a number of position/ velocity updates. 




Finally, we update the temporally smoothed observables as 
follows: 



Vj,t+i = Vj,t + a(Vj,t - Vj,t)St. (33) 

Here St is the time between successive x 2 M2M steps. All the 
differential equations here are ordinary differential equations 
of the form dyi(t)/dt = fi(t, j/i, • • • , un), and the yi,n's in 
our case are the particle weights Wi or time-smoothed ob- 
servables yi at t n - We integrate them using a simple Euler 
method yi, n +i = y%,n + h f(t n ,yi,n) with t n+1 = t n + h 
and time step h — St. We could replace the E uler method 
by th e second-order Runge-Kutta method (c/. iPress et al.l 
1992), which is more expensive and requires more memory. 
Since we are not interested in the details of how the weights 
converge, but only in the final converged system, a simple 
Euler method suffices for our purposes. We write e in equa- 
tion (J31J as e = e'e" with e" = 10/ vaaXi^{Kji Aj/a(Yj)}. 
Thus e" times the last term in equation (|15[) is of order 



unity and we choose 2e' < a to avoid excessive temporal 
smoothing. 

The NMAGIC (N-particle Made-to-measure 
AlGorithm minimizing Chi squared) correction rou- 
tine can be combined with a standard A^-body code 
including a potential solver and time integrator, or a 
fixed-potential routine and integrator when the target is to 
be modeled in a given gravitational potential. This last case 
is most similar to the Schwarzschild method. In most of the 
tests below, we use a fixed potential expanded in spherical 
harmonics. 

However, in test E we allow the potential to vary, as we 
evolve from one triaxial model to another. For advancing the 
particles we use a standard leap frog time integrator with 
fixed time-step. The time-step value chosen leads to fluctu- 
ations of energy and angular momentum with amplitudes 
5 x 10~ 6 and 2 x 10~ 5 around their initial values, without 
systematic drift, over 80 half mass dynamical times in the 
fixed potential case. 

For test E, which models a triaxial system, a simple 
spherical harmonic expansion suffices for s olving for the po - 
tential. We follow the method described bv lSellwoodl < |2003h : 
we tabulate coefficients of a spherical harmonic expansion of 
the density on a 1-D radial grid but retain the exact angular- 
dependence up to some adopted Z ma x, the maximum order 
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of the spherical harmonic expansion. We include terms up 
to i m ax = 4 in this experiment. Particles are binned on th e 
radial grid using the scheme described bv lSellwoodl (|2003h . 
This then gives the forces on the grid, from which we in- 
terpolate back to a particle's position for the gravitational 
forces. Test E involves a cuspy model; in order to properly 
resolve this we use a radial grid at radii = e 74 — 1 with 
7 = ln(r max + l)/£ 

in a x 

; we use £ 

max — 301 for 301 shells and 

r ma x = 40. 

NMAGIC is written in Fortran 90 and parallelized with 
the MPI library. We distribute the N particles as nearly 
evenly as possible over N p processors. Parallelizing in only 
the observables would not scale well with large N p , because 
of the different nature of the observables, and would require 
a large memory on each processor when N is large. In Figure 
Q] we present a high level flowchart of the operational logic 
of NMAGIC. 

In order to test the scaling with Np of NMAGIC we 
considered N = 1.8 X 10 6 and N = 816 observables (640 
density and 176 kinematic) with N p varying from 1 to 120. 
These values of N p and N axe adequate for the experiments 
presented here and are used in test C of Table[T] Since we are 
only interested in the scaling of the x 2 M2M parallelization 
with N p , we only execute the x 2 M2M 50 times, without 
recomputing the potential or advancing particles. In Figure 
[2]we present these scaling results as time per step (left hand 
axis, plus symbols) and steps per unit time (right hand axis, 
open squares) as functions of N p . We generally find that 
our implementation of x 2 M2M scales very well with N p . 
Defining the speedup S(N P , N) as 



S(N P ,N) 



T(1,N) 
T(N P ,N) 



(34) 



where T(N P , N) is the time for computing N particles on N v 
processors, we fit a standard Amdahl's law ( Amda hj |l967h 



S(N P ,N) 



f + (l-f)/N p ' 



(35) 



in order to determine the fraction of sequential code, /. We 
obtained that / ~ 0.010, i.e. the sequential part of the code 
accounts for only 1%. 



4 TARGET MODELS AND THEIR 
OBSERVABLES 

We will test the NMAGIC code on spherical, axisymmetric, 
and triaxial target models. The spherical target is a particle 
model constructed from the analytic density and distribu- 
tion function of an isotropic Hernquist sphere. As oblate 
target we take a maximally rotating three-integral model. 
Finally, we construct both a stationary and a rotating tri- 
axial target system. We use the NMAGIC code itself to gen- 
erate a dynamical equilibrium structure for these models. It 
will be seen that the x 2 M2M method provides a very useful 
means to set up dynamical equilibrium models of galaxies 
for which no analytic distribution functions are known, in 
order to study the properties of such systems. 

In the following subsections, we describe in turn each 
of these targets and their construction. We determine the 
target observables obtained from these models, and describe 
how we obtain errors for these observables. These will be 



50 



100 



Figure 2. The performance of our implementation of x 2 M2M. 
We used 1.8x 10 s particles without potential calculations or parti- 
cle motion. On the left hand axis we label time per step required, 
with the corresponding data indicated by plus symbols, while on 
the right hand axis we label steps per unit time, with the corre- 
sponding data now shown by open squares. Note that the scale 
is logarithmic on the left and linear on the right. The fraction of 
sequential code, /, from these data was computed at ~ 1%. 



needed in Section [S] where we present the results of building 
X 2 M2M models to match these targets. The reader who is 
mainly interested in these tests of NMAGIC can in a first 
reading directly go to that section. 



4.1 Spherical Target 

Our fi rst target is a spherical isotropic Hernquist (|Hernquistl 
Il990h model, which we will refer to as target SIH. Its density 
and potential are given by 

aM , , GM 



g(r) 



2nr(r + a) 3 



ip(r) 



r + a 



(36) 



where a is the scale length, M is the total mass, and G 
is the gravitational constant. The projected effective (half- 
mass) radius equals _R c h « 1.8153a. We use units such that 
M = a = 1. The target mass Mk on shell rk is given by the 
sum of the contributions of the adjacent shells, 



M k = 4tt / g(r)-yk IC (r)r 2 dr. 



(37) 



The innermost (outermost) shell is an exception because 
only the layer immediately exterior (interior) contributes. 

We construct SIH models on a radial grid with 40 shells, 
quasi-logarithmically spaced in radius with inner and outer 
boundaries at r m i n = 5 x 10~ 4 and r max = 20. The distri- 
bution function is truncated at -E max = 0(r max ). At that 
truncation, the mass included is 



Aftrun 



Bmi„ = V(0) 



dM 
~d~E 



dE, 



(38) 



with (dM/dE) the diffe rential energy distribution {e.g. 
iBinnev fc Tremaind [l987h and thus Mt run = 0.86. Figure 
[3] compares the mass on shells (hereafter "mass profile") 
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Figure 3. The mass in shells profile computed from equation J36D 
is shown as the solid line, whereas the dashed line illustrates the 
mass profile computed from a spherical Hernquist particle model 
generated from a truncated DF. 
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Figure 4. The upper panel shows the mass profile computed 
from equation (1406 for q = 0.6 (solid line) and from a Hernquist 
particle model made from a DF and squeezed along the z-axis 
(dashed line). The lower panel is the same but for A2Q. 



Mp(rk) for a particle realization of this truncated distribu- 
tion function (constructed u sing the method described in 
iDebattista fc Sellwoodl [2000 ) , with the Mk from the Hern- 
quist density profile as in equation (|37|l . For small radii 
the mass profiles match but for larger radii Mp(rk) is 
significantly smaller than Mk due to the finite extent of 
the particle realization, consisting only of particles with 
E < iJmax = ^(r max ). Using Mk as target observables would 
increase the mass of particles on the outer (near) circular 
orbits and would therefore increase the tangential velocity 
dispersion. We will thus use the Mp(rk) as targets and omit 
the subscript P in the following. We also include zero-valued 
higher order mass moments to enforce sphericity. 

We assume Poisson errors for the radial mass: a(Mj) = 
\J MjMtrun/N where N is the total number of particles used 
in the particle model. For the errors in the higher order 
mass moments, we use Monte-Carlo experiments in which 
we generate particle realizations of the density field of the 
target model using 1.8 x 10 6 particles, which is the same 
number as in the x 2 M2M models. 

Kinematics of the targe t can be comput e d from a DF. 
We u se the isotropic DF l|Hernquistl Il990l . ICarollo et al.l 
1 19951 ) 



f{E) 



■(3arcsin(g) + q(l — g 2 ) 1 ^ 2 



(l_ g 2)5/2 

x(l-2g 2 )(8g 4 -8g 2 -3)) 



(39) 



with g = y/-aE/GM, and E is the energy. We determine 
kinematic observables of the target on a projected radial 
grid with 30 shells, quasi-logarithmically spaced in radius 
and bounded by i? m i n = and i? ma x = 10 = 5.51i? c ff- On 
the shell midpoints we compute the and /14 moments of 
the isotropic Hernquist model from the DF of equation (1391) . 
We will use integral field-like kinematic data to recover the 



spherical targets in Section [5] More precisely, we multiply 
the /12,/s and /u,fc moments by the projected mass of the 
truncated SIH model within each radial grid shell to ob- 
tain the mass-weighted higher order moments Mk hz,h an d 
Mk /u,fc, which we use as the target observables. While this 
procedure is not perfectly self-consistent, because the mo- 
ments are from the infinite extent analytic DF while the 
mass is from the truncated DF, the differences are very 
small. The main advantage of doing this is that it allows 
us to compute the uncertainties in these kinematic observ- 
ables, which we assume a(Mkh nt k) = a(h n ) M c -J Mk /M c 
with a(hn) — 0.005, Mk the target mass in shell k, and M c 
the mass in the central grid shell. 

4.2 An Oblate Three-Integral Target made with 
NMAGIC 

Our oblate target model has density 
aM 



g(m) 



2-Kqm(ra + a) 3 



(40) 



where M and a are total mass and scale radius, and m 2 = 
R 2 + {z/q) 2 with q being the flattening. This density belongs 
to th e family of flattened 7 models (|Dehnen fc Gerhard! 
1 19941 ). wi th 7=1. We compute th e gravitational potential 
from (cf. iBinnev fc Tremainelll987l . section 2.3) 



<p(R,z) 



with 



GM 4>(m)dr 
2a Jo (l+r)^r + g 2 



R 2 



+ 



T + 1 ' T + q 2 ' 

m 2 + 2am 
(m + a) 2 



(41) 

(42) 
(43) 
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Figure 5. Kinematic major and minor axis slits for the oblate 
models, with the cells along each slit indicated. The ellipse cor- 
responds to one ij c ff of the equivalent Hcrnquist model squeezed 
by 0.6 in the z direction. 

by numerical integration, and tabulate it using a coarse and 
a fine linear grid in the meridional (R — z) plane. The coarse 
grid extends to R = z = 30a with 500 x 500 grid points. 
To increase the resolution at small 7? and z we replace the 
20x20 "innermost" grid cells at (R, z) = (0, 0) to (1.2a, 1.2a) 
by a finer grid also consisting of 500 x 500 grid points. 

In our experiments, we view the model edge-on along 
the x-axis as line-of-sight. Our targets are the mass moments 
Aim,k of the three-dimensional density p, and - for these 
oblate models - the kinematic moments m h n , n = 1, 4. 
We define an effective radius R c a m 1.8153a which is equal 
to that of the spherical Hernquist model. We set M = a = 1 
and q = 0.6. The target mass moments Ai mt k on shell rh are 
given by the sum of the contributions of the adjacent shells 
and are computed through 

Ai m>k = J Y lm (e,^)g(^ IC (r)d 3 x. (44) 

The innermost (outermost) shell is an exception because 
only the layer immediately exterior (interior) contributes. 
The setup of the radial grid is identical to that used for the 
spherical model and for our tests below we use Mj,, A2o,k, 

A22,k, ^66, fc- 

Figure U compares Mk and A%q : k computed from equa- 
tion (|40|) with Mp(ru) and Ap t 2o{rk) obtained from a spher- 
ical Hernquist particle realization built from a DF and 
squeezed along the z-axis by q = 0.6. As in Figure^ Mp(r k ) 
and Ap t 2o(rk) match Mk and A2o,k within r < 5a but then 
approach zero at larger radii towards r max . This difference 
is again due to the finite extent of the particle model. Be- 
low we therefore use the radial mass profile Mp and the 
higher mass moments Ap t i m as targets, and again we omit 
the subscript P in the following. 

We assume errors in the target mass profile a{Mj) as 
for the spherical model. For the errors in the higher order 



Figure 6. Target mass and A20 profiles for the triaxial models. 
The solid line shows target T54 while the dashed line shows target 
T53. 

mass moments, we use Monte-Carlo experiments in which 
particle realizations of the density field of the target model 
are generated using 5 x 10 5 particles, which is the same 
number as in the x 2 M2M models. 

In our oblate models we attempt to recover the target 
system from both slit and integral field kinematic data. Thus 
as kinematic target observables we use the projected mass- 
weighted Gauss-Hermite moments along the major and mi- 
nor axes in Test C, and on a grid of 30 x 20 points covering 
positions on the sky in [—3.6,3.6] x [—1.8, 1.8] in Test D. A 
schematic representation of the slit setup is shown in Figure 
\5\ The slits extend out to about 2R e g ~ 3.6. 

The target kinematics are determined from a 4 x 
10 6 particle representation of a maximally rotating three- 
integral model for the density distribution of equation (|40|l 
with q — 0.6. This is constructed by first evolving an 
isotropic spherical Hernquist model to the desired shape, 
using x 2 M2M, and then switching the in-plane velocity vec- 
tors of all particles with positive angular momentum J z 
to negative J z , leading to a DF which is still a vali d so- 
lution of the Boltzmann equation (|Lvnden-B cll I960]). For 
each slit or integral field cell p we obtain the mass in 
that cell M p and the mass-weighted Gauss-Hermite mo- 
ments M p hi lP , ■ ■ ■ , Mp fi4, p . We assume errors for the mass- 
weighted Gauss-Hermite coeffici ents as for the spherical 
model: a(M p h n , p ) — a(h n )M c \J M p /M c , where M p is the 
mass in slit cell p, computed by Monte-Carlo integration. In 
this case, we set a(h n ) = 0.005 (0.003) for the central slit 
(integral field grid) cell m c to approximate realistic errors. 

4.3 Making Triaxial Models with NMAGIC 

In the tests below we also explore triaxial Hernquist target 
models with stellar densities 

B( s ) = 7i 7 , (45) 

v ' 2ttx z s(s + a) 3 v ' 
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Figure 7. Line-of-sight velocity field of the rotating triaxial parti- 
cle model (RT54K) as seen in the inertial frame. The co-rotation 
radius is R CO i ~ WR c f[. The FoV extends from — -R c ff to i? e fj 
along each direction. Its lower edge is parallel to the major axis, 
the line-of-sight is parallel to the intermediate axis. Notice the 
counter-rotation near the center. 



w here M is the total m ass, a the scale radius, and s = 
yj (x/xq) 2 + y 2 + (z/zo) 2 . Here y is the longest axis, and 
the parameters xq and z§ are the axis ratios. As before, we 
use units with M = a = 1 and we define the effective ra- 
dius with reference to the spherical model, i.e. R c g sa 1.815. 
We generate two targets with different triaxialities, charac- 
terized by the tria xiality parameter T = (1 — Xq)/(1 — Zq) 
|Franx et al.lll99lT ). The more triaxial target, hereafter T53, 
has xo — 0.9 and zo — 0.8 (T = 0.53) whereas the less 
triaxial target, hereafter T54, has xo = 0.85 and zo = 0.7 
(T = 0.54). In both cases the target is observed along its 
intermediate (:r-)axis. 

Like our oblate target model, the triaxial models cannot 
be represented by a DF based on the integrals of motion. 
We therefore construct them through particle realizations 
via a two step process. Starting from a spherical Hernquist 
particle realization made from a DF as before, we squeeze 
this along the x- and z-axes by factors xo and zo, respec- 
tively, and compute the desired target density observables 
Mk and the higher order mass moments A20, fc, A22, k up to 
A&o,k using the same radial binning as in the spherical and 
oblate targets. Ai m components with I > 6 are small and 
we omit them. The squeezing is rigid, i. e. without regard to 
the internal motions. We repeat this 30 times, squeezing the 
spherical Hernquist model rigidly along random orientations 
to the desired shapes. From these 30 particle representations 
of the model we compute the means and one a variations 
around the mean for the Ai m ^. The former are taken as 
target density observables, the latter as their errors. The 
uncertainties on the radial mass in shells profile are taken 
to be a(M k ) — M k Mt T un/N as before. Figure [6] shows the 
target mass and A20 profiles as functions of radius for T54 



(solid line) and T53 (dashed line) as well as their uncertain- 
ties. 

After this first step, which only gives target density 
observables, we then use x 2 M2M to evolve a spherical 
Hernquist model to generate self-consistent triaxial parti- 
cle realizations of T54 and T53. In addition we generate 
a slowly tumbling version of T54 with corotation radius 
R cor w 10-Rcff, by applying x 2 M2M in the appropriately 
rotating frame. The final models now have self-consistent 
kinematics; in order to distinguish them from the purely 
density targets we refer to them as models T53K and T54K 
for the non-tumbling models and RT54K for the tumbling 
model. 

These final self-consistent models T54K and RT54K 
can now be used as targets in their own right, and we 
can compute (observer frame) target kinematics m p h n , P 
from them. We compute the kinematics of both T54K and 
RT54K on a 12 x 12 grid extending from —Reg to R e s- 
For the uncertainties in the kinem atic observables we adopt 
a(m p h„ yP ) = a{h n )M c ^/ M p /M c with a(h n ) = 0.005 the er- 
ror in h n , M p the mass in grid cell p, and M c the mass in 
the central grid cell. The M p 's were obtained directly from 
the particles. The velocity field of the target system RT54K 
in the observer's frame is shown in Figure [7] This veloc- 
ity field is characterized by disk-like counter-rotation close 
to the mid-plane and near cylindrical rotation away from 
the plane. These kinematics for this slowly tumbling triaxial 
model represent a valid dynamical model, but are unlikely 
to be the unique dynamical solution for the model's density 
distribution. 



5 TESTS OF NMAGIC 

In this section we will use the x 2 M2M algorithm to solve 
some modeling problems of increasing dimensionality and 
complexity, starting with spherical systems and ending with 
rotating triaxial models. The goal of these experiments is 
to investigate the convergence of the code, the quality with 
which various data are modeled, and the degree to which 
known properties of the target models can be recovered from 
their simulated data. We will see how these issues depend on 
the initial model, geometry, and amount of data available. 

Table Q] lists all the experiments that we have carried 
out, including the target and the initial model identifica- 
tions. We will refer to the final x 2 M2M models by the prefix 
'F' to the test model name (e.g., FA for the final model of 
Test A). Generally, these final models are obtained in two 
steps. First we use only the target density observables in 
the x 2 M2M algorithm, and once these have converged, we 
add the kinematic observables. Finally, we integrate all or- 
bits for some time in the potential without x 2 M2M correc- 
tions to test whether equilibrium has been reached. Unless 
mentioned otherwise, we use 1.8 x 10 6 particles and set the 
entropy parameter p to a small (<§C 1) value; see the dis- 
cussion in Section 5.1.1. In most experiments, the particle 
distribution is evolved in the fixed target potential (this is 
analogous to the Schwarzschild modelling approach) , but we 
include one test (model E) in which we also let the gravita- 
tional potential evolve. 
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Test 


ICS 


Target 


e' 


e" 


//. 

r* 


A 


RP 


SIH 


0.025 


6.32 10~ 7 


4.3 10" 4 


A2 


RP 


SIH 


0.025 


1.32 10~ 6 


4.3 10 2 


A3 


RP 


SIH 


0.050 


1.24 10" 6 


4.3 10" 4 


A4 


RP 


SIH 


0.100 


6.80 10~ 7 


4.3 10" 4 


B 


SIH-2 


SIH 


0.025 


1.76 10~ 6 


4.3 10" 4 


C 


ORIH 


031 


0.05 


3.94 10" 7 





D 


ORIH 


031 


0.05 


3.94 10~ 7 





E 


T53K 


T54K 


0.15 


5.06 lO" 8 


4.3 10 2 


F 


T54K 


RT54K 


0.15 


3.77 10-* 


4.3 10 2 



Table 1. Tests of NMAGIC carried out in this paper, with model 
names and parameters. For all models, we have used a = 2.1e'. 



5.1 Spherical Models 

5.1.1 Initial model and time-evolution 

The aim of our first experiment, Test A, is to reproduce a 
spherical isotropic Hernquist (SIH) model by a 1.8 x 10 B par- 
ticle model. We start by generating a P lummer model from 
its DF {e.g. Binnev fc Tremaind 11987!) . u sing the method 
described in iDebattista fc Sellwoodl (|200Ch . The DF of the 
Plummer model is truncated at <J?(r max ), with r max = 20, 
and has a scale length b — 1 and unit total mass. We then re- 
lax these particles in the analytic Hernquist potential, which 
is held fixed while the particle orbits are integrated. We re- 
fer to the resulting particle distribution as initial model RP 
(relaxed Plummer). 

Then with x 2 M2M we first adjust the density distribu- 
tion of model RP to that of the target SIH, using as target 
observables M k — v47nAoo h (equation [37)) and A\ mk = 
for 1 < Z ^ 6, < m ^ Z (equation I23|) with Monte Carlo 
errors estimated as described in Section 14.11 After conver- 
gence the even kinematic moment observables Mfc/l2,& and 
MfcZi4 j fc are added with errors given also in Section [4.11 Fi- 
nally the system is integrated for some time without apply- 
ing the x 2 M2M corrections. 

The second experiment B is identical to A except that 
instead of model RP we use a second Hernquist model SIH- 
2 as initial conditions for NMAGIC. SIH-2 differs from the 
target model SIH in that its radial scale length a = 1.4 
instead of a — 1. 

Figure |8^ shows the time evolution of x 1 /N a of the 
particle model A during and after the x 2 M2M evolution. 
Throughout N refers to all the observables, density and 
kinematics, regardless of whether they are being used in the 
FOC or not; thus N is a constant. The time evolution of 
a sample of 100 particle weights of the SIH particle model 
is presented in Figure [8p. From these figures one sees that 
the overall x 2 /N decreases quickly at the beginning of both 
phases (density adjustment only phase, and a density and 
kinematic observable adjustment phase). However, particle 
weights keep evolving for significantly longer time-scales. For 
this reason we integrate and adjust particle weights in both 
phases for relatively long times, about 15 dynamical times 
at r max ■ 
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Figure 8. (a) Top: Time evolution of x 2 in test A. (b) Bottom: 
Time evolution of a set of 100 particle weights in test A. wq is 
the initial weight of the particles; wq = l/N. The time-interval 
plotted includes a first phase of density adjustment (t ^ 2250), a 
second phase of density and kinematic adjustment (2250 < t $C 
4500), and a final phase of free evolution during which the weights 
do not change (t ^ 4500). Time is in units where the dynamical 
time at the half-mass radius is 6.0, and the dynamical time at 
r ma x is 150. 



5.1.2 Convergence to the target observables for different 
initial conditions 

The fit of the final particle models FA and FB to the ob- 
servables is illustrated in Figure [9] The top panel shows the 
radial mass and A20 coefficient, whereas the bottom panel 
shows the kinematic targets and final model observables for 
m h,2 in the upper and m /14 in the lower panel. As can also 
be seen from Fig. |5J the final model fits the input data to 
within la. The corresponding error bars are smaller than 
the crosses in the top panel of Fig. [9] see Fig.[6]for an exam- 
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Figure 9. (a) Top: Radial mass and A20 profiles for the target 
model SIH and the final models FA from Test A (solid line) and 
FB from Test B (dashed line), (b) Bottom: Kinematic profiles 
m/12 and m/14, for the same models. - In all panels, the data 
points with errors correspond to the SIH target, the solid line 
to the final particle model FA, and the dashed line to the final 
model FB. The error bars in the target mass distribution are not 
shown as they are smaller than the symbol sizes. The absolute 
errors shown decrease outwards due to the mass weighting; the 
corresponding relative errors increase outwards. 

pie. The same is true for A20 except when the target values 
are zero as in Fig. [9] Error bars for the mass observables are 
therefore not plotted in this and subsequent similar figures. 

All model observables in Fig.[9]are temporally smoothed 
observables as in equation (1 1 Q p . After some free evolution 
with x 2 M2M turned off both models fit the target data 
within the errors. The free evolution is necessary because 
X 2 M2M pushes the model towards a perfect fit to the ob- 
servables, at the expense of continually changing particle 
weights. Deviations are largest in the outer parts where or- 
bital time-scales are longest. Model FB, which had an ini- 
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-0.5 0.5 

Log(w/w ) 



Figure 10. Histogram of the particle weights in the final FA 
model, obtained from Plummer model initial conditions (solid 
line). The dashed line shows the histogram of particle weights 
when spherical Hernquist ICs with scale length a = 1.4 were used 
(FB). uiq is the initial weight of the particles; u>o = 1/JV in all 
cases. 



tial particle distribution closer to the target, is generally 
smoother and fits the data better, but differences are within 
the errors. NMAGIC achieved satisfactory models even from 
the less favorable, cored Plummer initial conditions. 

Figure \W\ compares the histogram of final particle 
weights of the FA and FB models, all normalized by their 
initial weight. Model FA has a significant tail towards high 
weights, and a peak at correspondingly lower particle weight 
such that the mean particle weight is the same as for the 
more symmetric weight distribution of model FB. On aver- 
age, the weights of particles in model FA had to change by 
more than those in model FB. We can quantify this by defin- 
ing an effective particle number Neg characterizing mass 
fluctuations through 

N cft =N^, (46) 
w 

where w and w 2 are the mean and mean-square particle 
weights. This reduces to N for equal-mass particles, to one 
when one particle dominates, and discards particles with 
near-zero weights. For the final models FA and FB the ef- 
fective numbers of particles are AT eff = 5.7xl0 5 and 1.5xl0 6 , 
respectively, while for both models N = 1.8 x 10 6 . 

The origin of this difference between the two models 
can be seen from Figure II lb . which plots the radial den- 
sity profile of the target SIH (stars), the initial models RP 
and SIH-2, and the temporally smoothed final models. We 
computed the densities using the identical radial grid as was 
used for the mass targets. The density profile of the SIH tar- 
get is well reproduced by the final particle models FA and 
FB across more than a factor of 100 in radius. The largest 
relative deviation in the density 5p/p occurs at small radii 
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Figure 11. (a) Top: Radial density profiles in the spherical mod- 
els. Uppermost panel: Density profiles for the Hernquist target 
profile, SIH (stars), the final models FA (dashed line) and FB 
(dash-dotted line), and their respective initial condition models 
RP and SIH-2 (dotted and dash-triple-dotted lines). Middle panel: 
Relative deviation from the target density Ap/p, for the two mod- 
els FA and FB using the same line styles, (b) Bottom: Differential 
energy distributions. The truncated analytic Hernquist DF used 
for target SIH is shown by the star symbols. The dashed line 
corresponds to the final x 2 M2M model FA, and the dotted line 
indicates the relaxed Plummer initial conditions RP. 



Figure 12. (a) (b) Top: Internal kinematics of the final model 
FA. The upper panel show ay, a v and erg, the lower panel the 
v r , V(p and vg . The stars correspond to the analytic oy from the 
untruncated DF. Model FA is very nearly isotropic and has negli- 
gible rotation, despite starting from anisotropic initial conditions. 
Bottom: Anisotropic internal kinematics of the initial model RP. 
The dotted, dashed, and dash-dotted lines show ay, (Tip, and ag 
of the RP particle model. For comparison, the solid line corre- 
sponds to the analytic oy of the untruncated analytic DF of the 
SIH target model. 



and never reaches more than 5%. In this region, model RP 
has few particles and the large relative error is due to Pois- 
son noise. Model FB, which starts out closer to the target 
SIH fits better in this region. 

Model RP is clearly significantly less dense than SIH 
inside r ~ 0.3a; it has a core whereas the target profile is 
cuspy. Also, it has a steeper outer density profile than the 
target model. To match model RP to SIH therefore requires 



NMAGIC to increase the particle masses both in the central 
regions and in the outer halo of the model. This causes the 
high- weight tail in the distribution in Fig. 1101 as we verified 
by inspecting the positions of particles with Wi > 2wo- 

Figure HTb presents the differential energy distributions. 
The final particle model FA matches the analytic differential 
energy distribution of the isotropic Hernquist model (equa- 
tion 1391) very well. 
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As a final test, Figure 112b shows the intrinsic veloci- 
ties (lower panel) and velocity dispersions (upper panel) of 
the analytic, untruncated DF and the final x 2 M2M model 
FA. The match to the target kinematics is good and model 
FA is nearly isotropic, despite the fact that it has evolved 
from an initial RP model that is moderately anisotropic. 
The anisotropy of the initial model RP is shown in Figure 
[T2b which compares its intrinsic velocity dispersions o>, o~ v 
and erg with the analytic ay of the SIH target model. The 
residual anisotropy in model FA is caused by the relative 
absence of radial orbits resulting from truncating the DF. 

5.1.3 Dependence on e and fi 

In the tests described so far, we have used e' = 0.025 for the 
correction steps in the FOC. In general, small values of e' 
result in a smooth evolution but slow convergence, whereas 
large values of e' change the global model too rapidly to at- 
tain a properly phase-mixed stationary solution. Thus gen- 
erally we have found e < 0.1 to give good results. This is 
illustrated in Fig. 1131 which shows that test A converges to 
essentially identical density distributions and differential en- 
ergy distributions for values of 0.025 ^ e' ^ 0.1 (models FA, 
FA3, FA4). Only for the largest value e' = 0.1 do we start 
seeing small deviations in the density profile of more than a 
few percent from the target model. Also, the effective parti- 
cle number [equation (|46|l ] decreases from 5.7 x 10 5 through 
3.3 x 10 5 to 1.0 x 10 5 for models FA, FA3, and FA4, respec- 
tively. Thus we will generally use e < 0.1, but because the 
speed of convergence also depends on the number and kind 
of observables used for the corrections, we have sometimes 
also increased e slightly. Figure [14] shows the distributions 
of particle weights for these models. They develop larger 
wings for larger values of e' . Because particles weights are 
then changed by larger amounts, the reshuffling is greater 
until convergence is reached. 

In models FA and FB, we have also set the entropy pa- 
rameter /jtoa small (<IC 1) value, which allows the NMAGIC 
code to concentrate on fitting the data. (Note that, because 
the term Kij Aj /a(Yj) in the FOC is large, even fi — 1.0 
leads to only a small contribution of the entropy terms in 
the FOC). While the purpose of not setting /j, to zero exactly 
originally was to prevent overly large fluctuations in the par- 
ticle weights, in fact, a test with /j, = has given essentially 
identical results to the ones reported. Fig. [13] shows that 
also for model FA2 with 10 6 times larger entropy parameter 
than in model FA, the target density and differential en- 
ergy distribution are fitted equally well as before. Generally, 
the best value to use for the entropy depends on the initial 
model, the data to be fitted, and the intrinsic structure of 
the target, and it must be determined separately for each 
application. A more systematic investigation of the effect of 
the entropy term is therefore deferred to a future paper in 
which we will use x 2 M2M to model and determine mass-to- 
light ratio, anisotropy, etc., for a real galaxy. 

5.2 Oblate Models 

The task we set the algorithm here is a difficult one: start- 
ing from a non-rotating system, we see whether we can re- 
cover the maximally rotating three-integral model described 
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Figure 13. (a) Top: Radial density profiles for various spherical 
models constructed for the Hernquist target profile, SIH. Upper 
panel: Density profiles for the target model (stars), the model 
FA (dashed line) and several tests that differ from model FA 
by the values of the parameters e' and u (see Table 1). Middle 
panel: Relative deviation from the target density Ap/p, for the 
same models, (b) Bottom: Differential energy distributions. Stars: 
target model SIH. Lines: same models as in top panel. 



in Section 14.21 in which the weights of all counter-rotating 
particles should be zero. We perform two such experiments, 
one using slit data as kinematic targets (Test C), the sec- 
ond using integral field kinematic targets (Test D). As in 
the spherical experiments, we keep the potential fixed while 
evolving the system with x 2 M2M in runs C and D. 

Both experiments start from an initial model which is 
constructed by relaxing a spherical Hernquist particle model 
consisting of 5 x 10 s particles in the oblate potential. As in 
experiments A and B, we then apply x 2 M2M in 2 steps, first 
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Figure 14. Histogram of the particle weights in the final FA 
model, obtained from Plummer model initial conditions (solid 
line). The other histograms show the particle weight distributions 
for models FA2 (dashed), FA3 (dash-dotted), and FA4 (dotted). 
wo is the initial weight of the particles; wq = 1/iV in all cases. 
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Figure 16. Mass- weighted higher order moments along the major 
and minor axes for the slit-reconstructed oblate model FC. The 
target observables are shown as error bars, whereas the model 
observables for model FC are indicated by the dashed lines, re- 
spectively. Kinematics along the major axis are shown on the left 
and those along the minor axis on the right. 
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Figure 15. The mass and mass A20 profiles for the oblate models. 
The data points show the target and the lines shows the converged 
models FC (dashed) and FD (full). 

for the density alone, and when this has converged, for both 
the density and kinematics. The density part of the runs is 
identical for experiments C and D. 

Figure [15] plots the mass and A20 radial profiles of the 
target (error bars) and the final x 2 M2M models FC and FD. 
As in the spherical tests, the target density distribution is 
very well fitted by the x 2 M2M models. 

The mass- weighted kinematics along the major and mi- 



nor axes of model FC are shown in Figure [T6l while Figure 
[T7]show the as-observed kinematics of both models. The lat- 
ter are calculated by dividing the mass-weighted moments 
by the mass in the slit resp. grid cell, and using the rela- 
t ions v = ftarg — V^ctargfei and a — a t!iIg — v / 2o" ta r g /i2 (e.g. 
, iRix et all Il997t ). All kinematic quantities for the recon- 
structed models are shown At — 500 (> 3 dynamical times 
at r max ) after switching off the x 2 M2M corrections. The 
fits are generally excellent except for the higher order mo- 
ments near the boundaries of the kinematic fit regions, where 
counter-rotating particles with high energies still make sig- 
nificant contributions, because their weights have not yet 
been sufficiently reduced. 

Figure [18] showing the weight distributions for both 
models FC and FD clearly illustrates the stronger con- 
straints placed on the model by the integral field data. 
In both models, the NMAGIC code works at reducing the 
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Figure 19. Particle weight distributions projected onto the (E, L 
the initial relaxed isotropic Hernquist model (top right), and the tw 
bottom right) and from density and integral field kinematics (bottc 

weights of the counter-rotating particles, but has clearly 
gone a lot further in model FD. 



) plane, for the maximally rotating three-integral target (top left), 
a models reconstructed from density and slit kinematic targets (FC, 
tn left). 



5.3 Triaxial Models 

5.3.1 Evolving the potential self- consistently 



Finally, in Figurc[T!|]we show the distribution of weights 
in the (E,L Z ) plane for the target model, initial relaxed 
model, and the two models FC and FD. The success of the 
X 2 M2M method in removing the counter-rotating particles 
amply present in the initial model is apparent, particularly 
for model FD. Of course, in applications aimed at obtaining 
a best-fit representation of some galaxy kinematic data it 
would have been smart to start the iterations from an initial 
model that is better adapted to the problem at hand. 



We illustrate NMAGIC 's capabilities with two very differ- 
ent triaxial model experiments. In run E, we start with 
the self-consistent model T53K as initial conditions and use 
NMAGIC to converge to target T54K. With this model, we 
test the full capabilities of x 2 M2M, which make this tech- 
nique more general than Schwarzschild's method: in model 
E, we solve for the potential as the system evolves and follow 
the model in its self-consistent potential throughout, akin to 
an TV-body experiment. For this purpose we use the spher- 
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Figure 17. Gaussian best fit velocity (top left), velocity disper- 
sion (top right), Gauss-Hermite moments /13 (bottom left) and /14 
(bottom right) along the major axis (left) and minor axis (right), 
for the models with slit data targets (dashed line), integral field 
kinematic targets (full), and the initial model (dotted). The error 
bars show the target kinematics. 
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Figure 18. Distribution of particle weights in the final models 
FC (top) and FD (bottom). 



ical harmonic potential solver described in Section [3] above 
and update the potential after every 25 x 2 M2M steps. 

The resulting final model FE gives an excellent match 
to the density of the target model T54K, as is apparent from 
comparing the Mk and A20 profiles in Figure 1201 Figure [2X1 
shows the kinematics within R c s of the models T54K and 
FE. All mass-weighted kinematic observables mhi,..., m/2,4 
of the final model match the target observables at better 
than one a over almost the entire FoV, except for a few iso- 
lated regions reaching two a. The random location of these 
deviations imply that they are due only to Poisson noise in 
the target model, the observables of which have not been 
temporally smoothed. 



5.3.2 Rotating vs. non-rotating models 

Test F is an interesting experiment in different ways. Start- 
ing from T54K, we use NMAGIC to attempt to converge to 
the observables of the tumbling target model RT54K, with 
a triaxial model which does not tumble but remains station- 
ary relative to the observer. Thus this experiment explores 
whether it is possible to identify a kinematic signature of 
slow figure rotation in elliptical galaxies. Since the initial 
conditions possess neither rotation nor internal net stellar 
streaming, if this model fails to converge it may well be be- 
cause the problem admits no solution. Because of this, test 
F is interesting in its own right, apart from as a validation 
of NMAGIC. 

In fact, NMAGIC was able to converge the mass- 
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Figure 20. Mass and A20 profiles for experiment E. The dots 
show the target T54K and the solid lines show FE. 




Figure 21. The difference between kinematics in model T54K 
and model FE. The observables of FE are the temporally 
smoothed mass-weighted moments while those of T54K are not 
temporally smoothed. The differences have been divided by the 
corresponding assumed errors. The FoV extends from — -R e ff to 
R c f{ along each direction. 



weighted kinematic moments to within about one a of their 
target values; however, the residuals maps (Figure [22]) show 
spatially correlated residuals in mh\ . When we compare the 
global velocity field of model FF with that of RT54K we 
find that the degree of cylindrical rotation around the tum- 
bling axis (z-axis) is higher in RT54K than it is in model 
FF (Figure I23[) . Near the mid-plane, instead, the velocity 
field of both models is very similar, including the counter- 
rotation seen near the center. We can explore whether the 




Figure 22. The difference between kinematics in target RT54K 
and model FF. The observables are the mass- weighted higher or- 
der moments, and have been divided by the corresponding as- 
sumed errors. The kinematics of RT54K are instantaneous but 
those of FF are time-averaged. The FoV extends from — i? e ff to 
R e ff along each direction. 



residual differences are due to having assumed too large er- 
rors in the mass- weighted moments by decreasing the errors 
by a factor of five. The corresponding final model looks very 
similar to model FF but now with reduced x 2 > 4. Thus the 
difference is likely intrinsic and can be used to recognize a 
tumbling galaxy. A more complete analysis of this problem 
will be undertaken elsewhere (De Lorenzi et al. in progress) . 



6 CONCLUSIONS 

We have presented a made-to-measure algorithm for con- 
structing particle models of stellar systems from observa- 
tional data, building o n the made-to-measure method of 
ISver fc Tremainel (|l99fj| . ST96) . An important element of our 
new method is the use of the standard x 2 merit function at 
the heart of the algorithm, in place of the relative error used 
by ST96. The improved algorithm, which we label x 2 M2M, 
allows us to assess the quality of a model for a set of tar- 
get data directly, using a statistically well-defined quantity 
(x 2 )- Moreover, this quantity is well-defined and finite also 
when a target observable takes on zero values. 

This property has enabled us to incorporate kinematic 
observables including higher-order Gauss-Hermite moments 
into the force-of-change equation. Kinematic and density (or 
surface density) observables can then be used simultane- 
ously to correct the particle weights. The price of changing 
to x 2 M2M from the original formulation is that the kernels 
which project the particle weights and phase-space coordi- 
nates into model observables cannot themselves depend on 
the particle weights. In general this is quite natural for (vol- 
ume or surface) density observables. For the kinematics this 
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Figure 23. Left: Line-of-sight velocity field of the final non- 
tumbling triaxial particle model FF. Right: Difference of the line- 
of-sight velocity fields between the non-tumbling triaxial parti- 
cle model FF and the tumbling target galaxy RT54K divided by 
the errors as described in the tex t. We a ssume an error in the 
mean velocity cr{vj) = \f2cr(h\ )^J m c /mj(Tj , where we assumed 

a {hi) = - yjl/2cr{v) / a l|Rix et al.lll997l ). In both panels the FoV 
extends from —R a g to i? e ff along each direction. 



means that we need to use mass- weighted kinematic observ- 
ables. Nonetheless, this is not a significant limitation. 

We have implemented the x 2 M2M method in a fast, 
parallel code, NMAGIC. This code also incorporates an op- 
tional but fast potential solver, allowing the potential to 
vary along with the model density. Its implementation of 
the x 2 M2M algorithm is highly efficient, with a sequential 
fraction of only ~ 1%. This has allowed us to build various 
models with large numbers of particles and based on many 
observables, and to run them for ~ 10 6 steps. 

Then we have carried out a number of tests to illustrate 
the capabilities and performance of NMAGIC, employing 
spherical, oblate and triaxial target models. The geometric 



flexibility by itself is one of the main strengths of the method 
- no symmetry assumptions need to be made. 

In the spherical experiments NMAGIC converged to 
the correct isotropic model from anisotropic initial condi- 
tions, demonstrating that a unique solution, if present, can 
be recovered. Both the truncated distribution function and 
the intrinsic velocity dispersions were recovered correctly. 
Two initial models with different density distributions were 
used in these experiments. While both converged to the final 
isotropic model, that with density closer to the density of 
the final model had smaller final deviations from the target 
observables, and a narrower distribution of weights. In both 
experiments, the observables (density and integral field-like 
kinematics) each converged in a few dynamical times at the 
outer boundary td,o, whereas the particle weights kept evolv- 
ing for significantly longer, ~ 10td,o- 

In the oblate experiments we gave the algorithm a dif- 
ficult problem to solve. The target system was a maximally 
rotating three-integral model in which the weights of all 
counter-rotating particles were zero. Using density observ- 
ables and either slit or integral field kinematics, NMAGIC 
was asked to recover this maximally rotating model start- 
ing from an isotropic spherical system relaxed in the oblate 
potential. After about lOO'OOO correction steps, particle 
weights on the counter-rotating side were reduced by a fac- 
tor of ~ 50, the distribution of weights approached that of 
the target, and a good fit to the kinematic constraint data 
was achieved. Only near the boundary of the kinematic data 
did particles on orbits further out, whose weights had not 
yet converged, still cause some deviations from the target 
kinematics. These experiments also clearly showed the ad- 
vantage of integral field data over slit data for constraining 
the model. 

Our triaxial experiments showed that it is possible to 
start from one triaxial model and converge to another. We 
anticipate that this ability will be very useful in constructing 
models for the triaxial elliptical galaxies with which nature 
confronts us. One of these triaxial experiments included a 
potential update step every 25 x 2 M2M steps, demonstrating 
that including an evolving potential is also practical. 

In the final experiment, we first generated a particle 
model of a slowly tumbling triaxial system to use as a tar- 
get. We then matched its volume density and line-of-sight 
kinematics with a stationary model. We showed that the 
mass- weighted kinematic moments of the figure rotating sys- 
tem was fitted to within one a by the non-rotating system 
out to -Reff- However the residuals in the first order kine- 
matic moment are correlated, which gives a clear signature 
of tumbling which the non-tumbling model is not able to 
match, even when the assumed errors are decreased by a 
significant factor. We thus conclude that, at least for this 
triaxial system, it is possible to distinguish between internal 
stellar streaming and pattern rotation within _R e ff provided 
a full velocity field is available. A more complete study of 
this problem will be presented elsewhere. 

This experiment also demonstrates the usefulness of the 
X 2 M2M algorithm for modeling mock (rather than real) 
galaxies in order to learn about their dynamics. We note 
that such an experiment would not have been practical with 
standard iV-body simulations. 

Compared to the Schwarzschild method, the main 
advantages of the x 2 M2M algorithm as implemented in 
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NMAGIC are that (i) stellar systems without symmetry 
restrictions can be handled relatively easily, (ii) it avoids 
complicated procedures for sampling, binning, and storing 
orbits, and (iii) the potential can be evolved self-consistently 
if needed. In the examples given, a simple isotropic spher- 
ical model was evolved into a suitable initial model, which 
contained the required wide range of orbital shapes. Every 
X 2 M2M model corresponds to a new set-up of a complete 
orbit library in the Schwarzschild method; so in problems 
where the same orbit library can be reused, Schwarzschild's 
method will be faster. However, NMAGIC is highly parallel, 
so suites of models with ~ 10 6 particles are feasible on a PC 
cluster. 

There is clearly room for improving the current imple- 
mentation of the x 2 M2M algorithm, and there is a need 
to study carefully the parameters that enter the algorithm, 
such as magnitude and frequency of the correction steps, 
entropy, etc., which we will address in future work. 

However, the different applications presented in this pa- 
per show that the x 2 M2M algorithm is practical, reliable and 
can be applied to various dynamically relaxed systems. High 
quality dynamical models of galaxies can be achieved which 
match targets to ~ la for plausible uncertainties in the ob- 
servables, and without symmetry restrictions. We conclude 
that x 2 M2M holds great promise for unraveling the nature 
of galaxies. 
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